import sys
sys.path.append('C:/Users/Hu/Dropbox/Research/PythonWork/Cancer/src/STAT/')

import numpy as np
import csv
from datetime import datetime
from time import strftime
import matplotlib.pyplot as plt
import math
import random
from scipy import stats
import ols
from scipy.stats import f


def getCurTime():
    """
    get current time
    Return value of the date string format(%Y-%m-%d %H:%M:%S)
    """
    format='%Y-%m-%d %H:%M:%S'
    sdate = None
    cdate = datetime.now()
    try:
        sdate = cdate.strftime(format)
    except:
        raise ValueError
    return sdate

def build_data_list(inputCSV):
    sKey = []
    fn = inputCSV
    f = open(inputCSV)
    #ra = csv.DictReader(file(fn), dialect="excel")
    ra = csv.DictReader(f, dialect="excel")
    
    for record in ra:
        #print record[ra.fieldnames[0]], type(record[ra.fieldnames[-1]])
        for item in ra.fieldnames:
            temp = int(float(record[item]))
            sKey.append(temp)
    sKey = np.array(sKey)
    sKey.shape=(-1,len(ra.fieldnames))
    return sKey

def geneBpChoices(tmin, tmax, num):
    tmin = float(tmin)
    tmax = float(tmax)
    num = int(num)
    step = (tmax - tmin)/num
    tbpchoices = []
    for i in range(num):
        tbpchoices.append(tmin+i*step)
    return tbpchoices

def pickBreakpoint(response, predictor, bpvarible, bpChoices, predictorName):
    results = np.zeros((len(bpChoices)-1, 2))
    print bpChoices
    bpvarible = np.array(bpvarible)
    predictor = np.array(predictor)
    predictor_t = np.transpose(predictor)
    predictorName = np.append(predictorName, 'bpvarible')

    for i in range(len(bpChoices)-1):
        x2star = (bpvarible - bpChoices[i]) * np.greater(bpvarible, bpChoices[i])
        tempPredictor = np.append(predictor_t, x2star)
        tempPredictor.shape = (len(predictorName),-1)
        tempPredictor = np.transpose(tempPredictor)
        #print tempPredictor
        tempmodel = ols.ols(response, tempPredictor,'y', predictorName)
        results[i,0] = i
        results[i,1] = tempmodel.R2
    #print results

    optBP = int(results[np.argmax(results, axis = 0)[1],0])
    print 'Optimal Index:', optBP
    print 'Optimal changepoint: ', bpChoices[optBP], ' exp value: ', np.exp(bpChoices[optBP]), ' with R2 = ', results[optBP, 1]

    x2star = (bpvarible - bpChoices[optBP]) * np.greater(bpvarible, bpChoices[optBP])
    tempPredictor = np.append(predictor_t, x2star)
    tempPredictor.shape = (len(predictorName),-1)
    optPredictor = np.transpose(tempPredictor)
    optmodel = ols.ols(response, optPredictor,'y', predictorName)
    #optmodel = ols.ols(response, optPredictor,'y',predictorName)
    
    y_hat = np.dot(optmodel.x, optmodel.b)
    #fileLoc = filepath + 'pieceR_model2_y_hat.csv'
    #np.savetxt(fileLoc, y_hat, delimiter=',', fmt = '%s')
    
    print optmodel.summary()
    print 'MSE =', optmodel.sse
    #print 'before bp, b0 =',optmodel.b[0], ', b1 =', optmodel.b[1], ', bd =', optmodel.b[2]
    #print 'after bp, b0 =', optmodel.b[0] - optmodel.b[3] * bpChoices[optBP], ', b1 =', optmodel.b[1], ', bd =', optmodel.b[2] + optmodel.b[3]

    #calpredictedvalue(zip(x1, predictor), bpChoices[optBP], zip(optmodel.b, [optmodel.b[0] - optmodel.b[3] * bpChoices[optBP], optmodel.b[1], optmodel.b[2] + optmodel.b[3]]), 'exp_inoutflow_model2A.csv')
    return results, y_hat    

def separateData(data1, data2, dataC, threshold):
    if (len(data1) <> len(data2)) or (len(data1) <> len(dataC)):
        print 'Error in input data length! --- seperateData.error.1'

    data1Left = []
    data1Right = []
    data2Left = []
    data2Right = []
    dataleftIdx = []
    datarightIdx = []
    
    for i in range(len(data1)):
        if dataC[i] < threshold:
            data1Left.append(data1[i])
            data2Left.append(data2[i])
            dataleftIdx.append(i)
        else:
            data1Right.append(data1[i])
            data2Right.append(data2[i])
            datarightIdx.append(i)

    return np.array(data1Left), np.array(data1Right), np.array(data2Left), np.array(data2Right), dataleftIdx, datarightIdx

def calculateR2(tleftmodel, trightmodel, totalmean):
    #self.sse = dot(self.e,self.e)/self.df_e
    #tsse = np.dot(tmodel.e, tmodel.e)/tmodel.df_e
    #tsse = np.dot(tmodel.e, tmodel.e)
    tsse = tleftmodel.sse * tleftmodel.df_e + trightmodel.sse * trightmodel.df_e
    
    #y_hat = np.dot(tmodel.x,tmodel.b)
    tsst = np.dot(tleftmodel.y - totalmean, tleftmodel.y - totalmean)
    tsst = tsst + np.dot(trightmodel.y - totalmean, trightmodel.y - totalmean)
    
    #print tleftmodel.sse * tleftmodel.df_e, trightmodel.sse * trightmodel.df_e, np.dot(tleftmodel.y - totalmean, tleftmodel.y - totalmean), np.dot(trightmodel.y - totalmean, trightmodel.y - totalmean)
    #print tleftmodel.y.var() * tleftmodel.df_e, trightmodel.y.var() * trightmodel.df_e
    #print 1 - (tleftmodel.sse * tleftmodel.df_e)/np.dot(tleftmodel.y - totalmean, tleftmodel.y - totalmean)
    #print 1 - (trightmodel.sse * trightmodel.df_e)/np.dot(trightmodel.y - totalmean, trightmodel.y - totalmean)
    
    return 1-tsse/tsst
    #print ((tleftmodel.nobs+trightmodel.nobs-1.0), (tleftmodel.nobs+trightmodel.nobs-tleftmodel.ncoef-trightmodel.ncoef)), ((tleftmodel.nobs+trightmodel.nobs-1.0)/(tleftmodel.nobs+trightmodel.nobs-tleftmodel.ncoef-trightmodel.ncoef))
    #return 1 - (tsse/tsst)*((tleftmodel.nobs+trightmodel.nobs-1.0)/(tleftmodel.nobs+trightmodel.nobs-tleftmodel.ncoef-trightmodel.ncoef))

def simpleRegression(response, predictor, predictorName):
    model = ols.ols(response, predictor,'y',predictorName)
    #model = ols.ols(response, predictor)
    print model.summary()
    print model.b
    print predictor
    #print 'b0 =',model.b[0], ', b1 =', model.b[1], ', bd =', model.b[2], ', MSE =', model.sse
    return np.dot(model.x, model.b)

def separateRegression(response, predictor, sepData, bpChoices, predictorName, equaltestid):
    results = np.zeros((len(bpChoices)-1, 4))
    print bpChoices
    chosencriterion = 'r2'  # max(r2) is the same criterion as min(sse)
    for bpid in range(len(bpChoices)-1):
        #print bpid
        responseLeft, responseRight, predictorLeft, predictorRight, dataleftIdx, datarightIdx = separateData(response, predictor, sepData, bpChoices[bpid])
        #print responseLeft, responseRight, predictorLeft, predictorRight
        leftmodel = ols.ols(responseLeft, predictorLeft,'y',predictorName)
        rightmodel = ols.ols(responseRight, predictorRight,'y',predictorName)
        results[bpid,0] = bpid
        if chosencriterion == 'r2':           
            results[bpid,1] = leftmodel.R2adj
            results[bpid,2] = rightmodel.R2adj
            #results[bpid,3] = 1 - (leftmodel.e.var() + rightmodel.e.var())/(leftmodel.y.var() + rightmodel.y.var())
            results[bpid,3] = calculateR2(leftmodel, rightmodel, np.mean(response))
            #results[bpid,3] = (leftmodel.R2 + rightmodel.R2)/2
        elif chosencriterion == 'sse':
            results[bpid,1] = leftmodel.sse * leftmodel.df_e
            results[bpid,2] = rightmodel.sse * rightmodel.df_e
            results[bpid,3] = results[bpid,1] + results[bpid,2]
            #yhatL = np.dot(leftmodel.x, leftmodel.b)
            #print results[bpid,1], np.sum((responseLeft - yhatL)**2)

    #print results
    if chosencriterion == 'r2':
        optBP = int(results[np.argmax(results, axis = 0)[-1],0])
    elif chosencriterion == 'sse':
        optBP = int(results[np.argmin(results, axis = 0)[-1],0])

    responseLeft, responseRight, predictorLeft, predictorRight, dataleftIdx, datarightIdx = separateData(response, predictor, sepData, bpChoices[optBP])
    leftmodel = ols.ols(responseLeft, predictorLeft,'y',predictorName)
    rightmodel = ols.ols(responseRight, predictorRight,'y',predictorName)
    #print responseLeft, predictorLeft,'y',predictorName

    #equaltestid = int(equaltestid)
    if equaltestid[0] > -1:
        temppredictorLeft = predictorLeft
        temppredictorLeft[:,equaltestid[1]] = temppredictorLeft[:,equaltestid[0]] + temppredictorLeft[:,equaltestid[1]]
        temppredictorLeft = np.delete(temppredictorLeft, equaltestid[0], 1)
        temppredictorRight = predictorRight
        temppredictorRight[:,equaltestid[1]] = temppredictorRight[:,equaltestid[0]] + temppredictorRight[:,equaltestid[1]]
        temppredictorRight = np.delete(temppredictorRight, equaltestid[0], 1)
        temppredictorName = np.delete(predictorName, equaltestid[0], None)
        templeftmodel = ols.ols(responseLeft, temppredictorLeft,'y',temppredictorName)
        temprightmodel = ols.ols(responseRight, temppredictorRight,'y',temppredictorName)
        fleft = (leftmodel.R2 - templeftmodel.R2) * (leftmodel.nobs - len(predictorName) - 2) / (1 - leftmodel.R2)
        fright = (rightmodel.R2 - temprightmodel.R2) * (rightmodel.nobs - len(predictorName) - 2) / (1 - rightmodel.R2)
        pleft = 1 - f.cdf(fleft, 1, leftmodel.nobs - len(predictorName) - 2)
        pright = 1 - f.cdf(fright, 1, rightmodel.nobs - len(predictorName) - 2)

    yhatL = np.dot(leftmodel.x, leftmodel.b)
    yhatR = np.dot(rightmodel.x, rightmodel.b)
    yhat = np.zeros(len(response))

    for i in range(len(yhatL)):
        yhat[dataleftIdx[i]] = yhatL[i]

    for i in range(len(yhatR)):
        yhat[datarightIdx[i]] = yhatR[i]

    yhat = np.exp(yhat)
    #fileLoc = filepath + 'separateR_model2_y_hat.csv'
    #np.savetxt(fileLoc, yhat, delimiter=',', fmt = '%s')
    print 'Optimal Index:', optBP
    print 'Optimal changepoint: ', bpChoices[optBP], ' exp value: ', np.exp(bpChoices[optBP]), ' with R2 = ', calculateR2(leftmodel, rightmodel, np.mean(response))
   
    print '----------------------------- left model -----------------------------'
    print leftmodel.summary()
    print '----------------------------- right model -----------------------------'
    print rightmodel.summary()

    print 'Optimal Index:', optBP
    print 'Optimal changepoint: ', bpChoices[optBP], ' exp value: ', np.exp(bpChoices[optBP]), ' with R2 = ', results[optBP, -1]
   
    outputstring = 'before bp'
    for i in range(len(predictorName)+1):
        outputstring += ', b' + str(i) + ' = ' + "%.2f" %(leftmodel.b[i]) + '(' + "%.2f" %(leftmodel.se[i])+ ')'
    outputstring += ', with R2 = ' + "%.2f" %(leftmodel.R2)
    if equaltestid[0] > -1:
        outputstring += ', p1 <> p2 with pvalue = ' + "%.2f" %(pleft)
    #print 'before bp, b0 =',leftmodel.b[0], ', b1 =', leftmodel.b[1], ', bd =', leftmodel.b[2], ', MSE =', leftmodel.sse
    print outputstring

    outputstring = 'after bp'
    for i in range(len(predictorName)+1):
        outputstring += ', b' + str(i) + ' = ' + "%.2f" %(rightmodel.b[i]) + '(' + "%.2f" %(rightmodel.se[i])+ ')'
    outputstring += ', with R2 = ' + "%.2f" %(rightmodel.R2)
    if equaltestid[0] > -1:
        outputstring += ', p1 <> p2 with pvalue = ' + "%.2f" %(pright)
    #print 'after bp, b0 =', rightmodel.b[0], ', b1 =', rightmodel.b[1], ', bd =', rightmodel.b[2], ', MSE =', rightmodel.sse
    print outputstring

    #calpredictedvalue(predictor, bpChoices[optBP], zip(leftmodel.b, rightmodel.b), 'exp_inoutflow_model2B.csv')
    #calconfidenceinterval(predictorLeft, predictorRight, [leftmodel.sse, rightmodel.sse], response, predictor, bpChoices[optBP], zip(leftmodel.b, rightmodel.b), 'ci_model2B.csv')
    return results, yhat

def calconfidenceinterval(leftx, rightx, mse, y, x, bp, b,filename):
    leftx = np.array(zip(np.ones(len(leftx)), leftx[:,0], leftx[:,1]))
    rightx = np.array(zip(np.ones(len(rightx)), rightx[:,0], rightx[:,1]))

    templeftx = np.linalg.inv(np.dot(leftx.T, leftx))
    temprightx = np.linalg.inv(np.dot(rightx.T, rightx))

    b = np.array(b)
    t = [stats.t.ppf(0.975, len(leftx) - len(leftx[0,:])), stats.t.ppf(0.975, len(rightx) - len(rightx[0,:]))]

    yhat = []
    i = 0
    #print len(x)
    for item in x:
        if item[1] < bp:
            bpflag = 0
            tempx = templeftx
        else:
            bpflag = 1
            tempx = temprightx
        temp = b[0, bpflag] + b[1, bpflag] * item[0] + b[2, bpflag] * item[1]
        xh = np.array([1, item[0], item[1]])
        #print np.dot(np.dot(xh, tempx), xh.T), mse[bpflag]*(np.dot(np.dot(xh, tempx), xh.T)
        se = math.pow(mse[bpflag]*(np.dot(np.dot(xh, tempx), xh.T)), 0.5) * t[bpflag]
        pise = math.pow(mse[bpflag]*(1 + np.dot(np.dot(xh, tempx), xh.T)), 0.5) * t[bpflag]
        cil = 0
        cih = 0
        pil = 0
        pih = 0
        if y[i] < (temp-se):
            cil = 1
        if y[i] > (temp+se):
            cih = 1
        if y[i] < (temp-pise):
            pil = 1
        if y[i] > (temp+pise):
            pih = 1
        yhat.append([temp, temp-se, temp+se, temp-pise, temp+pise,cil,cih,pil,pih])
        i += 1

    yhat = np.exp(yhat)
    #fileLoc = filepath + filename   #'exp_inoutflow_model2A.csv'
    headerstr = 'yhat,cileft,ciright,pileft,piright,cil,cih,pil,pih'
    #print fileLoc
    #np.savetxt(fileLoc, yhat, delimiter=',', header = headerstr, fmt = '%f')
    
    #plt.plot(range(len(y)), y, 'g', label='Observed')
    #plt.plot(range(len(y)), yhat[:,0], 'b', label='Predicted')
    #plt.plot(range(len(y)), yhat[:,1], 'b', label='CI Left')
    #plt.plot(range(len(y)), yhat[:,2], 'b', label='CI Right')
    #plt.plot(range(len(y)), yhat[:,3], 'r', label='PI Left')
    #plt.plot(range(len(y)), yhat[:,4], 'r', label='PI Right')
    #plt.show()

def cal_dist(x1, y1, x2, y2):
    temp_dist = (float(x1 - x2))**2 + (float((y1 - y2))**2)
    #print x1, y1, x2, y2, (float(x1 - x2))**2, (float((y1 - y2))**2), temp_dist
    temp_dist = temp_dist**0.5
    #print temp_dist

    return temp_dist

def calpredictedvalue(x, bp, b, filename):
    # x: [f1f2, dist]
    # b = [[b00, b10], [b01, b11], [b02, b12]]
    print b
    yhat = []
    b = np.array(b)
    bpflag = 1  # 0: before, 1: after
    for item in x:
        if item[1] < bp:
            bpflag = 0
        else:
            bpflag = 1
        temp = b[0, bpflag] + b[1, bpflag] * item[0] + b[2, bpflag] * item[1]
        #temp = np.dot(item, b[:,bpflag])
        yhat.append(temp)

    yhat = np.array(yhat)
    yhat = np.exp(yhat)
    #print yhat
    #fileLoc = filepath + 'model2b_yhat.csv'
    #np.savetxt(fileLoc, yhat, delimiter=',', fmt = '%f')
    expinoutflow = np.zeros((len(inoutflow), 2)) # [inflow, outflow, obs-exp inflow, obs-exp outflow]
    #print len(yhat), len(data)
    #print zip(data, yhat)
    i = 0
    for item in data:
        expinoutflow[int(item[1]), 1] += yhat[i]
        expinoutflow[int(item[2]), 0] += yhat[i]
        i += 1

    fileLoc = filepath + filename   #'exp_inoutflow_model2A.csv'
    headerstr = 'inflow,outflow'
    #np.savetxt(fileLoc, expinoutflow, delimiter=',', header = headerstr, fmt = '%f')
    
def cleardata(tempdata):
    #delete the row with nan or negative vlue
    #print len(tempdata)
    deletelist = []
    for i in range(len(tempdata)):
        #for j in range(len(tempdata[i])):
        j = -1
        if tempdata[i,j] < 0.0000001 or np.isnan(tempdata[i,j]):
            #print len(tempdata), i
            deletelist.append(i)
            #tempdata = np.delete(tempdata, i, 0)
    tempdata = np.delete(tempdata, deletelist, 0)
    #print tempdata
    return tempdata

def expbin(centroidCSV, grossCSV, bininterval, exp_dist_bin, graphCSV, exp_graphdist_bin):
    centroid = build_data_list(centroidCSV)
    grossflow = build_data_list(grossCSV)
    graphdist = build_data_list(graphCSV)
    totalflow = sum(grossflow[:,0])
    for i in range(len(centroid)):
        for j in range(len(centroid)):
            #if i <> j:
            if i > -1:
                dist = int(cal_dist(centroid[i, 0], centroid[i, 1], centroid[j, 0], centroid[j, 1]))
                temp = int(dist/bininterval)
                if temp > len(exp_dist_bin) - 1:
                    temp = len(exp_dist_bin) - 1
                exp = (float(grossflow[i,1]) * float(grossflow[j,0])) / totalflow
                exp_dist_bin[temp] += exp
                temp = graphdist[i,j]
                if temp > len(exp_graphdist_bin) - 1:
                    temp = len(exp_graphdist_bin) - 1
                exp_graphdist_bin[temp] += exp
    return exp_dist_bin, exp_graphdist_bin

def average_dist_based_on_graphdist(centroidCSV, graphCSV):
    centroid = build_data_list(centroidCSV)
    graphdist = build_data_list(graphCSV)
    output = np.zeros((100,2))
    for i in range(len(centroid)):
        for j in range(i+1,len(centroid)):
            tempdist = int(cal_dist(centroid[i, 0], centroid[i, 1], centroid[j, 0], centroid[j, 1]))
            tempgraphdist = graphdist[i,j]
            tempsum = output[tempgraphdist,0] * output[tempgraphdist,1] + tempdist
            output[tempgraphdist,0] += 1
            output[tempgraphdist,1] = tempsum/output[tempgraphdist,0]
    return output

def regression_with_bp(response, predictor, predictorName, bpvariate, bp):
    # bp: [breakpoint value (lnvalue), coe before bp, coe after bp]
    # all the data are lndata
    # response and bpvariate are one-dimensional
   # print response

    tsst = np.dot(response - np.mean(response), response - np.mean(response))
    
    for i in range(len(response)):
        if bpvariate[i] < bp[0]:
            bpcoe = bp[1]
        else:
            bpcoe = bp[2]

        response[i] = response[i] - bpvariate[i] * bpcoe

    #print np.var(response)*len(response)
    #print np.dot(response - np.mean(response), response - np.mean(response))

        
    responseLeft, responseRight, predictorLeft, predictorRight, dataleftIdx, datarightIdx = \
                  separateData(response, predictor, bpvariate, bp[0])
    #print dataleftIdx[:100], datarightIdx[:100]

    #print np.dot(responseLeft - np.mean(responseLeft), responseLeft - np.mean(responseLeft))
    #print np.dot(responseRight - np.mean(responseRight), responseRight - np.mean(responseRight))
    
    leftmodel = ols.ols(responseLeft, predictorLeft,'y',predictorName)
    rightmodel = ols.ols(responseRight, predictorRight,'y',predictorName)

    yhatL = np.dot(leftmodel.x, leftmodel.b)
    yhatR = np.dot(rightmodel.x, rightmodel.b)
    yhat = np.zeros(len(response))

    for i in range(len(yhatL)):
        yhat[dataleftIdx[i]] = yhatL[i]

    for i in range(len(yhatR)):
        yhat[datarightIdx[i]] = yhatR[i]

    #yhat = np.exp(yhat)
    print '----------------------------- left model -----------------------------'
    print leftmodel.summary()
    print '----------------------------- right model -----------------------------'
    print rightmodel.summary()

    tsse = leftmodel.sse * leftmodel.df_e + rightmodel.sse * rightmodel.df_e
    
    print 'Whole R2 = ', calculateR2(leftmodel, rightmodel, np.mean(response)), 'or ',1-tsse/tsst, '(', tsse, tsst, ')'
    #print np.sum(response), len(response), np.sum(responseLeft), len(responseLeft), np.sum(responseRight), len(responseRight)
    #print np.mean(response), np.mean(responseLeft), np.mean(responseRight)
    for i in range(len(yhat)):
        if bpvariate[i] < bp[0]:
            bpcoe = bp[1]
        else:
            bpcoe = bp[2]
        yhat[i] = yhat[i] + bpvariate[i] * bpcoe

    #tempyhat = simpleRegression(response, predictor, predictorName)

    '''
    for f, b in zip(tempyhat, bpvariate):
        if b < bp[0]:
            bpcoe = bp[1]
        else:
            bpcoe = bp[2]
        f = f + b * bpcoe
    '''
    return yhat

#--------------------------------------------------------------------------
#MAIN

if __name__ == "__main__":
    print '===================================================='
    print "begin at " + getCurTime()
    
    filepath2009 = 'C:/_DATA/migration_census_2009/'
    file = filepath2009 + 'regressiondata.csv'
    #data2009 = build_data_list(file)
    
    filepath2000 = 'C:/_DATA/migration_census_2000/'
    file = filepath2000 + 'regressiondata.csv'    
    data2000 = build_data_list(file)    #[dist,graphdist,grossout,grossin,popout,popin,vol]
    #print data2000
    '''
    #print max(data[:,0]), min(data[:,0]), max(data[:,1]), min(data[:,1])
    #exp2000 = []
    #total2000 = sum(data2000[:,-1])
    #for item in data2000:
        #exp2000.append(float(item[2])*float(item[3])/total2000)

    #exp2009 = []
    #total2009 = sum(data2009[:,-1])
    #for item in data2009:
        #exp2009.append(float(item[2])*float(item[3])/total2009)

    bininterval = 10000 # 10km
    #print bininterval
    dist_bin2000 = np.zeros(int(max(data2000[:,0])/bininterval)+1)
    graphdist_bin2000 = np.zeros(int(max(data2000[:,1]))+1)

    exp_dist_bin2000 = np.zeros(int(max(data2000[:,0])/bininterval)+1)
    exp_graphdist_bin2000 = np.zeros(int(max(data2000[:,1]))+1)
    exp_dist_bin2000, exp_graphdist_bin2000 = expbin(filepath2000+'census_county_centroid.csv',
                              filepath2000+'grossflow.csv', bininterval,
                              exp_dist_bin2000, filepath2000+'cnty_graphdistance.csv',
                              exp_graphdist_bin2000)
    
    exp_dist_bin2009 = np.zeros(int(max(data2009[:,0])/bininterval)+1)
    exp_graphdist_bin2009 = np.zeros(int(max(data2009[:,1]))+1)
    exp_dist_bin2009, exp_graphdist_bin2009 = expbin(filepath2009+'census_county_centroid.csv',
                              filepath2009+'grossflow.csv', bininterval,
                              exp_dist_bin2009, filepath2009+'cnty_graphdistance.csv',
                              exp_graphdist_bin2009)
    #print a

    for item in data2000:
        temp = int(item[0]/bininterval)
        dist_bin2000[temp] += item[-1]
        temp = int(item[1])-1
        graphdist_bin2000[temp] += item[-1]

    avg_exp_dist_bin2000 = np.mean(exp_dist_bin2000)
    avg_exp_graphdist_bin2000 = np.mean(exp_graphdist_bin2000)

    dist_bin2009 = np.zeros(int(max(data2009[:,0])/bininterval)+1)
    graphdist_bin2009 = np.zeros(int(max(data2009[:,1]))+1)

    for item in data2009:
        temp = int(item[0]/bininterval)
        dist_bin2009[temp] += item[-1]
        temp = int(item[1])-1
        graphdist_bin2009[temp] += item[-1]
        
    avg_exp_dist_bin2009 = np.mean(exp_dist_bin2009)
    avg_exp_graphdist_bin2009 = np.mean(exp_graphdist_bin2009)

    adj_dist_bin2000 = dist_bin2000/(exp_dist_bin2000/avg_exp_dist_bin2000)
    adj_graphdist_bin2000 = graphdist_bin2000/(exp_graphdist_bin2000/avg_exp_graphdist_bin2000)

    adj_dist_bin2009 = dist_bin2009/(exp_dist_bin2009/avg_exp_dist_bin2009)
    adj_graphdist_bin2009 = graphdist_bin2009/(exp_graphdist_bin2009/avg_exp_graphdist_bin2009)

    tickfontsize = 10

    #plot histogram of dist.
    dist_bin2000 = np.array(zip(np.array(range(len(dist_bin2000)))*(bininterval/1000), dist_bin2000))
    adj_dist_bin2000 = np.array(zip(np.array(range(len(adj_dist_bin2000)))*(bininterval/1000), adj_dist_bin2000))
    
    dist_bin2009 = np.array(zip(np.array(range(len(dist_bin2009)))*(bininterval/1000), dist_bin2009))
    adj_dist_bin2009 = np.array(zip(np.array(range(len(adj_dist_bin2009)))*(bininterval/1000), adj_dist_bin2009))
    
    graphdist_bin2000 = np.array(zip(np.array(range(len(graphdist_bin2000)))+1, graphdist_bin2000))
    adj_graphdist_bin2000 = np.array(zip(np.array(range(len(adj_graphdist_bin2000)))+1, adj_graphdist_bin2000))

    graphdist_bin2009 = np.array(zip(np.array(range(len(graphdist_bin2009)))+1, graphdist_bin2009))
    adj_graphdist_bin2009 = np.array(zip(np.array(range(len(adj_graphdist_bin2009)))+1, adj_graphdist_bin2009))
    
    dist_bin2000 = cleardata(dist_bin2000)
    adj_dist_bin2000 = cleardata(adj_dist_bin2000)
    dist_bin2009 = cleardata(dist_bin2009)
    adj_dist_bin2009 = cleardata(adj_dist_bin2009)
    graphdist_bin2000 = cleardata(graphdist_bin2000)
    adj_graphdist_bin2000 = cleardata(adj_graphdist_bin2000)
    graphdist_bin2009 = cleardata(graphdist_bin2009)
    adj_graphdist_bin2009 = cleardata(adj_graphdist_bin2009)   
    '''
    #inputPredictorName = ['lndist']
    #bpChoices = np.log(dist_bin2000[4:-2,0])
    #inputdata = adj_dist_bin2009
    
    #bpChoices = np.log(inputdata[4:-2,0])

    #reg = ols.ols(np.log(inputdata[1:,1]), np.log(inputdata[1:,0]),'y',inputPredictorName)
    #reg.summary()
    #result, yhat = separateRegression(np.log(inputdata[1:,1]), np.log(inputdata[1:,0]), np.log(inputdata[1:,0]),
                                      #bpChoices, inputPredictorName, [-2,1])
    #results, y_hat = pickBreakpoint(np.log(inputdata[1:,1]), np.log(inputdata[1:,0]), np.log(inputdata[1:,0]),
                                      #bpChoices, inputPredictorName)
    #disthist2000 = average_dist_based_on_graphdist(filepath2000+'census_county_centroid.csv',
                                                   #filepath2000+'cnty_graphdistance.csv')

    #disthist2009 = average_dist_based_on_graphdist(filepath2009+'census_county_centroid.csv',
                                                   #filepath2009+'cnty_graphdistance.csv')
    #headerstr = 'bin,cflow'
    #np.savetxt(filepath2000+'adj_graphdist_bin2000.csv', adj_graphdist_bin2000, delimiter=',', header = headerstr, fmt = '%s')

    data = data2000 #[dist,graphdist,grossout,grossin,popout,popin,vol]
    bp = [np.log(1000000), -1.335532, -0.374431]    # model9: adj_dist_2000
    #bp = [np.log(16), -1.686176, -0.732022]    # model10: adj_graphdist_2000
    #bp = [np.log(1040000), -1.432048, -0.306935]    # model9: adj_dist_2009
    #bp = [np.log(19), -1.780665, -0.498193]    # model10: adj_graphdist_2009

    lndata = np.log(data)
    print np.dot(lndata[:,-1] - np.mean(lndata[:,-1]), lndata[:,-1] - np.mean(lndata[:,-1]))

    #inputPredictor = []
    #for item in lndata:
        #inputPredictor.append([item[2], item[3], item[2]*item[2], item[3]*item[3]])
        #inputPredictor.append([item[2], item[3]])
    #inputPredictor = np.array(inputPredictor)
    inputPredictor = lndata[:,[2,3]]
    predictorName = ['lnf1','lnf2']
    #predictorName = ['lnf1', 'lnf2', 'lnf12', 'lnf22']
    #predictorName = ['lnp1','lnp2']
    yhat = regression_with_bp(lndata[:,-1], inputPredictor, predictorName, lndata[:,0], bp)
    #yhat = regression_with_bp(lndata[:,-1], lndata[:,[2,3]], predictorName, lndata[:,0], bp)
    #yhat = regression_with_bp(lndata[:,-1], lndata[:,[4,5]], predictorName, lndata[:,1], bp)
    yhat = np.exp(yhat)
    #print yhat
    
    headerstr = 'expflow'
    np.savetxt(filepath2000+'temp.csv', yhat, delimiter=',', header = headerstr, fmt = '%s')

    '''
    plotdata1 = dist_bin2000
    plotdata2 = dist_bin2009
    plotdata3 = graphdist_bin2000
    plotdata4 = graphdist_bin2009
    label1 = 'Original Data 2000'
    label2 = 'Original Data 2009'
    label3 = 'Adjusted Data 2000'
    label4 = 'Adjusted Data 2009'

    ax1 = plt.subplot(221)
    #plt.plot(dist_bin[:,0], dist_bin[:,1], marker = '^', linestyle = '', label = label1, ms=4, color ='w')
    #plt.plot(adj_dist[:,0], adj_dist[:,1], marker = 'o', linestyle = '', label = label2, ms=4, color ='w')
    plt.loglog(dist_bin2000[:,0], dist_bin2000[:,1], marker = '^', linestyle = '', label = label1, ms=4,
               markeredgewidth=1, markeredgecolor='b', markerfacecolor='None')
    plt.loglog(dist_bin2009[:,0], dist_bin2009[:,1], marker = 'o', linestyle = '', label = label2, ms=4,
               markeredgewidth=1, markeredgecolor='g', markerfacecolor='None')
    plt.xlabel('Distance (km)')
    plt.ylabel('Gross Flow')
    plt.xticks(fontsize = tickfontsize)
    plt.yticks(fontsize = tickfontsize)
    plt.legend(loc=3)
    #plt.axis([min(data[:83,0]/1000) * 0.8, max(data[:83,0]/1000) *1.5, min(data[:83,1]) * 0.8, max(data[:83,1]) * 3.5])
    
    ax2 = plt.subplot(222)
    #plt.plot(graphdist_bin[:,0], graphdist_bin[:,1], marker = '^', linestyle = '', label = label1, ms=4, color ='b')
    #plt.plot(adj_graphdist_bin[:,0], adj_graphdist_bin[:,1], marker = 'o', linestyle = '', label = label2, ms=4, color ='r')
    plt.loglog(graphdist_bin2000[:,0], graphdist_bin2000[:,1], marker = '^', linestyle = '', label = label1, ms=4,
               markeredgewidth=1, markeredgecolor='b', markerfacecolor='None')
    plt.loglog(graphdist_bin2009[:,0], graphdist_bin2009[:,1], marker = 'o', linestyle = '', label = label2, ms=4,
               markeredgewidth=1, markeredgecolor='g', markerfacecolor='None')
    plt.xlabel('Graph Distance')
    plt.ylabel('Gross Flow')
    plt.xticks(fontsize = tickfontsize)
    plt.yticks(fontsize = tickfontsize)
    plt.legend(loc=3)

    ax3 = plt.subplot(223)
    plt.loglog(adj_dist_bin2000[:,0], adj_dist_bin2000[:,1], marker = '^', linestyle = '', label = label3, ms=4,
               markeredgewidth=1, markeredgecolor='b', markerfacecolor='None')
    plt.loglog(adj_dist_bin2009[:,0], adj_dist_bin2009[:,1], marker = 'o', linestyle = '', label = label4, ms=4,
               markeredgewidth=1, markeredgecolor='g', markerfacecolor='None')
    plt.xlabel('Distance (km)')
    plt.ylabel('Gross Flow')
    plt.xticks(fontsize = tickfontsize)
    plt.yticks(fontsize = tickfontsize)
    plt.legend(loc=3)
    
    ax4 = plt.subplot(224)    
    plt.loglog(adj_graphdist_bin2000[:,0], adj_graphdist_bin2000[:,1], marker = '^', linestyle = '', label = label3, ms=4,
               markeredgewidth=1, markeredgecolor='b', markerfacecolor='None')
    plt.loglog(adj_graphdist_bin2009[:,0], adj_graphdist_bin2009[:,1], marker = 'o', linestyle = '', label = label4, ms=4,
               markeredgewidth=1, markeredgecolor='g', markerfacecolor='None')
    plt.xlabel('Graph Distance')
    plt.ylabel('Gross Flow')
    plt.xticks(fontsize = tickfontsize)
    plt.yticks(fontsize = tickfontsize)
    plt.legend(loc=3)
    
    plt.show()
        
    #lndata = np.log(data)

    #dist = data[:,0]
    #lndist = lndata[:,0]
    #graphdist = data[:,1]
    #lngraphdist = lndata[:,1]
    #bpChoices = np.log(range(min(graphdist)+2, max(graphdist)+1))
    #bpChoices = np.log(range(600000,4000000,100000))
    #bpChoices = np.log(range(20000,200000,20000))
    #bpChoices = np.log(range(600000,900000,100000))
    #bpChoices = np.log(range(40000,1000000,10000))
    #print bpChoices
    
    #inputPredictor = lndata[:,[2,3,1]]
    #inputPredictorName = ['logp1p2', 'logdist','oe','de']
    #inputPredictorName = ['lnp1', 'lnp2', 'lnp12', 'lnp22', 'lndist']
    #inputPredictorName = ['lnp1', 'lnp2', 'lnp12', 'lndist']
    #yhat = simpleRegression(output[:,0], inputPredictor, inputPredictorName)
    inputPredictorName = ['lnp1', 'lnp2', 'lndist']
    #result, yhat = separateRegression(lndata[:,-1], inputPredictor, lngraphdist, bpChoices, inputPredictorName, [0,1])

    inputPredictorName = ['lnp1', 'lnp2', 'lndist']
    #results, y_hat = pickBreakpoint(lndata[:,-1], inputPredictor, lndist, bpChoices, inputPredictorName)

    #fileLoc = filepath + 'piecewiseR_model2_results.csv'
    #fileLoc = filepath + 'separateR_model2_results.csv'
    headerstr = 'exp'
    #fileLoc = filepath + 'temp.csv'
    #np.savetxt(fileLoc, yhat, delimiter=',', header = headerstr, fmt = '%s')
    '''